Skip to content

Commit 40775ac

Browse files
santiagobadiasantiagobadia
authored
Merge pull request #801 from gridap/SpaceTimeByExtrusion
First steps for space time by extrusion (models with FE geo maps)
2 parents e27c6bb + 7aba09c commit 40775ac

14 files changed

+463
-1
lines changed

‎NEWS.md

+3
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
2121
## [0.17.13] - 2022-05-31
2222

2323
### Added
24+
- `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)
25+
- `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)
26+
- `GridWithFEMap`, which has a geometrical map defined by a FE function. Since PR [#801](https://github.com/gridap/Gridap.jl/pull/801)
2427
- Vertex bisection algorithm for refinement of triangular meshes in 2D. Since PR [#733](https://github.com/gridap/Gridap.jl/pull/733)
2528
- Generalized-α method for 1st order ODEs. Since PR [#781](https://github.com/gridap/Gridap.jl/pull/781)
2629
- Implemented (generalised) ModalC0 Polynomial bases and reference FEs. Since PR [#777](https://github.com/gridap/Gridap.jl/pull/777)

‎src/Arrays/Arrays.jl

+4
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,8 @@ export PosNegPartition
7878

7979
export FilterMap
8080

81+
export KeyToValMap
82+
8183
export MulAddMap
8284

8385
export Table
@@ -155,6 +157,8 @@ include("PosNegReindex.jl")
155157

156158
include("Reindex.jl")
157159

160+
include("KeyToValMaps.jl")
161+
158162
include("FilteredArrays.jl")
159163

160164
include("SubVectors.jl")

‎src/Arrays/KeyToValMaps.jl

+19
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
struct KeyToValMap{T<:Function} <: Map
2+
key_to_val::T
3+
end
4+
5+
function return_cache(m::KeyToValMap,key)
6+
K = typeof(key)
7+
V = typeof(m.key_to_val(key))
8+
d = Dict{K,V}()
9+
end
10+
11+
function evaluate!(cache,m::KeyToValMap,key)
12+
if haskey(cache,key)
13+
return get(cache,key,nothing)
14+
else
15+
val = m.key_to_val(key)
16+
cache[key] = val
17+
return val
18+
end
19+
end

‎src/FESpaces/AffineFEOperators.jl

+1-1
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
"""
33
AffineFEOperator
44
5-
Reprepresent a fully assembled affine (linear) finite element problem.
5+
Represent a fully assembled affine (linear) finite element problem.
66
See also [FEOperator](@ref)
77
"""
88
struct AffineFEOperator <: FEOperator

‎src/FESpaces/DiscreteModelWithFEMaps.jl

+149
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,149 @@
1+
"""
2+
Given a Discrete Model and a reffe, builds a new grid in which the geometrical
3+
map is a `FEFunction`. This is useful when considering geometrical maps that are
4+
the result of a FE problem (mesh displacement).
5+
"""
6+
struct GridWithFEMap{Dc,Dp,A,B,C,D,E} <: Grid{Dc,Dp}
7+
grid::Grid{Dc,Dp}
8+
fe_sp::A
9+
scal_fe_sp::B
10+
fe_map::C
11+
node_coords::D
12+
reffes::E
13+
end
14+
15+
function GridWithFEMap(model,order; kwargs...)
16+
orders = Fill(order,num_cells(model))
17+
# _reffes = [reffe]
18+
# reffes = lazy_map(Reindex(_reffes),cell_type)
19+
GridWithFEMap(model,orders; kwargs...)
20+
end
21+
22+
function GridWithFEMap(model,orders::AbstractArray; kwargs...)
23+
24+
# Create a FESpace for the geometrical description
25+
# The FEFunction that describes the coordinate field
26+
# or displacement can be a interpolation or a solution
27+
# of a mesh displacement problem
28+
T = eltype(get_node_coordinates(model))
29+
Ts = eltype(T)
30+
31+
os = orders
32+
_ps = get_polytopes(model)
33+
ct = get_cell_type(model)
34+
ps = lazy_map(Reindex(_ps), ct)
35+
36+
f(a,b) = LagrangianRefFE(T,a,b)
37+
reffes = lazy_map(f,ps,os)
38+
Vₕ = FESpace(model,reffes;conformity=:H1,kwargs...)
39+
40+
fs(a,b) = LagrangianRefFE(Ts,a,b)
41+
s_reffes = lazy_map(f,ps,os)
42+
Vₕ_scal = FESpace(model,s_reffes;conformity=:H1)
43+
44+
grid = get_grid(model)
45+
geo_map = get_cell_map(grid)
46+
47+
cell_ctype = get_cell_type(grid)
48+
c_reffes = get_reffes(grid)
49+
50+
# Create a fe_map using the cell_map that can be evaluated at the
51+
# vertices of the fe space (nodal type)
52+
# This returns a FEFunction initialised with the coordinates
53+
# But this is to build the FEFunction that will be inserted, it is
54+
# an advanced constructor, not needed at this stage
55+
56+
c_dofs = get_fe_dof_basis(Vₕ)
57+
dof_basis = get_data(c_dofs)
58+
c_nodes = lazy_map(get_nodes,get_data(c_dofs))
59+
60+
xh = zero(Vₕ)
61+
c_dofv = lazy_map(evaluate,dof_basis,geo_map)
62+
63+
Uₕ = TrialFESpace(Vₕ)
64+
65+
fv = get_free_dof_values(xh)
66+
dv = get_dirichlet_dof_values(get_fe_space(xh))
67+
gather_free_and_dirichlet_values!(fv,dv,Uₕ,c_dofv)
68+
69+
c_xh = lazy_map(evaluate,get_data(xh),c_nodes)
70+
c_scal_ids = get_cell_dof_ids(Vₕ_scal)
71+
72+
73+
nodes_coords = Vector{eltype(eltype(c_xh))}(undef,num_free_dofs(Vₕ_scal))
74+
Geometry._cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
75+
76+
GridWithFEMap(grid, Vₕ, Vₕ_scal, xh, nodes_coords, reffes)
77+
end
78+
79+
function _compute_node_coordinates(grid,xh)
80+
c_dofs = get_fe_dof_basis(grid.fe_sp)
81+
c_nodes = lazy_map(get_nodes,get_data(c_dofs))
82+
83+
c_xh = lazy_map(evaluate,get_data(xh),c_nodes)
84+
c_scal_ids = get_cell_dof_ids(grid.scal_fe_sp)
85+
86+
nodes_coords = grid.node_coords
87+
Geometry._cell_vector_to_dof_vector!(nodes_coords,c_scal_ids,c_xh)
88+
89+
return nothing
90+
91+
end
92+
93+
function add_mesh_displacement!(grid::GridWithFEMap,dh::FEFunction)
94+
95+
Xh = grid.fe_map
96+
97+
Xh_fv = get_free_dof_values(Xh)
98+
Xh_dv = get_dirichlet_dof_values(get_fe_space(Xh))
99+
100+
Xh_fv .= Xh_fv .+ get_free_dof_values(dh)
101+
Xh_dv .= Xh_dv .+ get_dirichlet_dof_values(get_fe_space(dh))
102+
103+
_compute_node_coordinates(grid,Xh)
104+
105+
end
106+
107+
function update_coordinates!(grid::GridWithFEMap,dh::FEFunction)
108+
109+
Xh = grid.fe_map
110+
111+
Xh_fv = get_free_dof_values(Xh)
112+
Xh_dv = get_dirichlet_dof_values(get_fe_space(Xh))
113+
114+
Xh_fv .= get_free_dof_values(dh)
115+
Xh_dv .= get_dirichlet_dof_values(get_fe_space(dh))
116+
117+
_compute_node_coordinates(grid,Xh)
118+
119+
end
120+
121+
# Interface
122+
get_node_coordinates(grid::GridWithFEMap) = grid.node_coords
123+
get_cell_node_ids(grid::GridWithFEMap) = get_cell_dof_ids(grid.scal_fe_sp)
124+
get_reffes(grid::GridWithFEMap) = grid.reffes
125+
get_cell_type(grid::GridWithFEMap) = get_cell_type(grid.grid)
126+
OrientationStyle(grid::GridWithFEMap) = OrientationStyle(grid.grid)
127+
RegularityStyle(grid::GridWithFEMap) = RegularityStyle(grid.grid)
128+
# santiagobadia: To think (does it make sense here, I think it is for surface problems)
129+
get_facet_normal(grid::GridWithFEMap) = get_facet_normal(grid.grid)
130+
131+
132+
# struct DiscreteModelWithFEMap{Dc,Dp} <: DiscreteModel{Dc,Dp}
133+
# model::DiscreteModel{Dc,Dp}
134+
# mapped_grid::GridWithFEMap{Dc,Dp}
135+
# function DiscreteModelWithFEMap(model::DiscreteModel{Dc,Dp},phys_map) where {Dc,Dp}
136+
# mapped_grid = GridWithFEMap(get_grid(model),phys_map)
137+
# new{Dc,Dp}(model,mapped_grid)
138+
# end
139+
# end
140+
141+
function DiscreteModelWithFEMap(model::DiscreteModel,order; kwargs...)
142+
mapped_grid = GridWithFEMap(model,order;kwargs...)
143+
MappedDiscreteModel(model,mapped_grid)
144+
end
145+
146+
function DiscreteModelWithFEMap(model,orders::AbstractArray; kwargs...)
147+
mapped_grid = GridWithFEMap(model,orders;kwargs...)
148+
MappedDiscreteModel(model,mapped_grid)
149+
end

‎src/FESpaces/FESpaces.jl

+16
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,15 @@ import Gridap.Arrays: autodiff_array_jacobian
4141
import Gridap.Arrays: autodiff_array_hessian
4242
import Gridap.Geometry: get_triangulation
4343
import Gridap.Geometry: get_cell_shapefuns
44+
import Gridap.Geometry: get_cell_type
45+
import Gridap.Geometry: MappedGrid
46+
import Gridap.Geometry: MappedDiscreteModel
47+
import Gridap.Geometry: get_node_coordinates
48+
import Gridap.Geometry: get_reffes
49+
import Gridap.Geometry: get_cell_node_ids
50+
import Gridap.Geometry: OrientationStyle
51+
import Gridap.Geometry: RegularityStyle
52+
import Gridap.Geometry: get_facet_normal
4453
import Gridap.CellData: attach_constraints_rows
4554
import Gridap.CellData: attach_constraints_cols
4655
import Gridap.CellData: CellField
@@ -194,6 +203,11 @@ export FESpaceWithLinearConstraints
194203

195204
export FiniteElements
196205

206+
export DiscreteModelWithFEMap
207+
export GridWithFEMap
208+
export add_mesh_displacement!
209+
export update_coordinates!
210+
197211
include("FESpaceInterface.jl")
198212

199213
include("SingleFieldFESpaces.jl")
@@ -240,6 +254,8 @@ include("DirichletFESpaces.jl")
240254

241255
include("FESpacesWithLinearConstraints.jl")
242256

257+
include("DiscreteModelWithFEMaps.jl")
258+
243259
export get_free_values
244260
function get_free_values(args...)
245261
@unreachable "get_free_values has been removed. Use get_free_dof_values instead."

‎src/Geometry/Geometry.jl

+4
Original file line numberDiff line numberDiff line change
@@ -124,6 +124,7 @@ export compute_reference_grid
124124

125125
export GridPortion
126126
export UnstructuredGrid
127+
export MappedGrid
127128

128129
export CartesianGrid
129130
export CartesianDescriptor
@@ -164,6 +165,7 @@ export CartesianDiscreteModel
164165
export GenericTriangulation
165166
export BoundaryTriangulation
166167
export DiscreteModelPortion
168+
export MappedDiscreteModel
167169

168170
export Interior
169171
export Boundary
@@ -214,6 +216,8 @@ include("UnstructuredDiscreteModels.jl")
214216

215217
include("CartesianDiscreteModels.jl")
216218

219+
include("MappedDiscreteModels.jl")
220+
217221
include("GridPortions.jl")
218222

219223
include("DiscreteModelPortions.jl")

‎src/Geometry/MappedDiscreteModels.jl

+81
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
function _cell_vector_to_dof_vector!(dof_vector,cell_node_ids, cell_vector)
2+
cache_cell_node_ids = array_cache(cell_node_ids)
3+
cache_cell_vector = array_cache(cell_vector)
4+
for k=1:length(cell_node_ids)
5+
current_node_ids = getindex!(cache_cell_node_ids,cell_node_ids,k)
6+
current_values = getindex!(cache_cell_vector,cell_vector,k)
7+
for (i,id) in enumerate(current_node_ids)
8+
dof_vector[current_node_ids[i]]=current_values[i]
9+
end
10+
end
11+
end
12+
13+
# """
14+
# MappedGrid
15+
16+
# Represent a grid with a geometrical map that is the composition of
17+
# a reference to a physical space map (standard)
18+
# and a (vector-valued) map from physical space to physical space. E.g.,
19+
# it can be used to include a high order map implemented by any map that is
20+
# a `CellField`.
21+
# """
22+
struct MappedGrid{Dc,Dp,T,M,L} <: Grid{Dc,Dp}
23+
grid::Grid{Dc,Dp}
24+
geo_map::T # Composition of old map and new one
25+
phys_map::M # New map in the physical space
26+
node_coords::L
27+
function MappedGrid(grid::Grid{Dc,Dp},phys_map) where {Dc,Dp}
28+
29+
@assert length(phys_map) == num_cells(grid)
30+
@assert eltype(phys_map) <: Field
31+
32+
function _compute_node_coordinates(grid,phys_map)
33+
cell_node_ids = get_cell_node_ids(grid)
34+
old_nodes = get_node_coordinates(grid)
35+
node_coordinates = Vector{eltype(old_nodes)}(undef,length(old_nodes))
36+
c_coor = get_cell_coordinates(grid)
37+
map_c_coor = lazy_map(evaluate,phys_map,c_coor)
38+
_cell_vector_to_dof_vector!(node_coordinates,cell_node_ids,map_c_coor)
39+
return node_coordinates
40+
end
41+
42+
model_map=get_cell_map(grid)
43+
geo_map=lazy_map(∘,phys_map,model_map)
44+
node_coords = collect(_compute_node_coordinates(grid,phys_map))
45+
new{Dc,Dp,typeof(geo_map),typeof(phys_map),typeof(node_coords)}(grid,geo_map,phys_map,node_coords)
46+
end
47+
end
48+
49+
function MappedGrid(grid::Grid{Dc,Dp},phys_map::Function) where {Dc,Dp}
50+
c_map = Fill(GenericField(phys_map),num_cells(grid))
51+
MappedGrid(grid,c_map)
52+
end
53+
54+
get_node_coordinates(grid::MappedGrid) = grid.node_coords
55+
get_cell_node_ids(grid::MappedGrid) = get_cell_node_ids(grid.grid)
56+
get_reffes(grid::MappedGrid) = get_reffes(grid.grid)
57+
get_cell_type(grid::MappedGrid) = get_cell_type(grid.grid)
58+
59+
"""
60+
MappedDiscreteModel
61+
62+
Represent a model with a `MappedGrid` grid.
63+
See also [`MappedGrid`](@ref).
64+
"""
65+
struct MappedDiscreteModel{Dc,Dp} <: DiscreteModel{Dc,Dp}
66+
model::DiscreteModel{Dc,Dp}
67+
mapped_grid::MappedGrid{Dc,Dp}
68+
function MappedDiscreteModel(model::DiscreteModel{Dc,Dp},phys_map) where {Dc,Dp}
69+
mapped_grid = MappedGrid(get_grid(model),phys_map)
70+
new{Dc,Dp}(model,mapped_grid)
71+
end
72+
end
73+
74+
get_grid(model::MappedDiscreteModel) = model.mapped_grid
75+
get_cell_map(model::MappedDiscreteModel) = get_cell_map(model.mapped_grid)
76+
get_grid_topology(model::MappedDiscreteModel) = get_grid_topology(model.model)
77+
get_face_labeling(model::MappedDiscreteModel) = get_face_labeling(model.model)
78+
79+
function Grid(::Type{ReferenceFE{d}},model::MappedDiscreteModel) where {d}
80+
get_grid(model)
81+
end

‎test/ArraysTests/KeyToValMapsTests.jl

+48
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
module KeyToValMapsTests
2+
3+
using Test
4+
using Gridap.Arrays
5+
using Gridap.ReferenceFEs
6+
7+
using Gridap
8+
using FillArrays
9+
10+
keys = [1,2,3,4,5]
11+
key_to_val(i) = i*ones(10)
12+
13+
m = KeyToValMap(key_to_val)
14+
r = lazy_map(m,keys)
15+
16+
r[1]
17+
c = return_cache(r)
18+
19+
_r = key_to_val.(keys)
20+
@test all(_r .== collect(r))
21+
22+
23+
key_to_val(i) = LagrangianRefFE(Float64,HEX,i)
24+
25+
m = KeyToValMap(key_to_val)
26+
keys = [1,2,3,4,5,4,3,2,1]
27+
r = lazy_map(m,keys)
28+
_r = key_to_val.(keys)
29+
@test all(get_order.(_r)-get_order.(r) .== 0)
30+
31+
end #module
32+
33+
# @time key_to_val(1)
34+
# @time for i in r end
35+
# @time for i in keys key_to_val(i) end
36+
37+
# @time key_to_val(5)
38+
# dict = Dict(keys .=> key_to_val.(keys))
39+
# # fetching from dict is fast
40+
# @time dict[5]
41+
# @time r[5]
42+
# a = getkey(dict,5,1)
43+
# a
44+
# return_cache(m,keys)
45+
46+
# haskey(dict,5)
47+
# get(dict,5,nothing)
48+
# get!(dict,5,nothing)

0 commit comments

Comments
 (0)