Grid refinement

Project.toml

@@ -13,6 +13,7 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"

src/Gridap.jl

@@ -56,6 +56,8 @@ include("MultiField/MultiField.jl")
 
 include("ODEs/ODEs.jl")
 
+include("Refinement/Refinement.jl")
+
 include("Exports.jl")
 
 end # module

src/Refinement/RefinedDiscreteModels.jl

@@ -0,0 +1,116 @@
+
+
+"""
+Refinement glue between two nested triangulations
+- `f2c_faces_map`     : Given a fine face gid, returns the gid of the coarse face containing it. 
+- `fcell_to_child_id` : Given a fine cell gid, returns the local child id within the coarse cell containing it.
+- `f2c_ref_cell_map`  : Ref coordinate map between the fcells and their parent ccell, i.e 
+                          for each fine cell defines Φ st. x_c = Φ(x_f)
+"""
+struct RefinementGlue{A,B,C} <: GridapType
+  f2c_faces_map::A
+  fcell_to_child_id::B
+  f2c_ref_cell_map::C
+
+  function RefinementGlue(f2c_faces_map,
+                          fcell_to_child_id,
+                          f2c_ref_cell_map)
+    A = typeof(f2c_faces_map)
+    B = typeof(fcell_to_child_id)
+    C = typeof(f2c_ref_cell_map)
+    new{A,B,C}(f2c_faces_map,fcell_to_child_id,f2c_ref_cell_map)
+  end
+end
+
+"""
+Ref coordinate map from fcell ref coords to ccell ref coords.
+Size -> Number of children per coarse cell.
+"""
+function get_f2c_ref_cell_map(reffe,ref)
+  ref_grid = Geometry.UnstructuredGrid(Visualization.compute_reference_grid(reffe,ref))
+  return get_cell_map(ref_grid)
+end
+
+"""
+Ref coordinate map between fine and coarse cells. 
+Size -> Number of fine cells. 
+"""
+function get_f2c_ref_coordinate_map(g::RefinementGlue)
+  m = Reindex(g.f2c_ref_cell_map)
+  return lazy_map(m,g.fcell_to_child_id)
+end
+
+#################################################################
+
+"""
+  RefinedDiscreteModel
+
+  `DiscreteModel` created by refining another `DiscreteModel`. 
+  The hierarchy is stored, allowing for the transfer of dofs 
+  between `FESpaces` defined on this model and its parent using a 
+  `RefinementTransferoperator`. 
+"""
+struct RefinedDiscreteModel{Dc,Dp,A,B} <: DiscreteModel{Dc,Dp}
+  model  ::A
+  parent ::B
+  glue   ::RefinementGlue
+
+  function RefinedDiscreteModel(model::DiscreteModel{Dc,Dp},parent::DiscreteModel{Dc,Dp},glue) where {Dc,Dp}
+    A = typeof(model)
+    B = typeof(parent)
+    return new{Dc,Dp,A,B}(model,parent,glue)
+  end
+end
+
+# DiscreteModel API
+Geometry.get_grid(model::RefinedDiscreteModel)          = get_grid(model.model)
+Geometry.get_grid_topology(model::RefinedDiscreteModel) = get_grid_topology(model.model)
+Geometry.get_face_labeling(model::RefinedDiscreteModel) = get_face_labeling(model.model)
+
+# Other getters
+get_model(model::RefinedDiscreteModel)  = model.model
+get_parent(model::RefinedDiscreteModel) = model.parent
+get_glue(model::RefinedDiscreteModel)   = model.glue
+
+function refine(model::DiscreteModel) :: RefinedDiscreteModel
+  @abstractmethod
+end
+
+
+# Cartesian builder 
+
+function RefinedCartesianDiscreteModel(domain::Tuple,nC::Int,ref::Int)
+  nF = ref*nC
+  # Models
+  parent = CartesianDiscreteModel(domain,(nC,nC))
+  child  = CartesianDiscreteModel(domain,(nF,nF))
+
+  # Glue 
+  faces_map      = [Int[],Int[],_create_f2c_cell_map(nC,ref)]
+  fcell_child_id = _create_child_map(nC,ref)
+  reffe          = LagrangianRefFE(Float64,QUAD,1)
+  ref_cell_map   = get_f2c_ref_cell_map(reffe,ref)
+  glue = RefinementGlue(faces_map,fcell_child_id,ref_cell_map)
+
+  # RefinedModel
+  model = RefinedDiscreteModel(child,parent,glue)
+  return model
+end
+
+function _create_f2c_cell_map(nC::Int,ref::Int)
+  nF = nC*ref
+
+  idx = Tuple.(CartesianIndices((nF,nF)))
+  a = map((i,j)->(1+(i-1)÷ref,1+(j-1)÷ref),first.(idx),last.(idx))
+  b = map((i,j)->(i-1)*nC+j,first.(a),last.(a))
+  return Array(reshape(transpose(b),nF*nF))
+end
+
+function _create_child_map(nC::Int,ref::Int)
+  nF = nC*ref
+  elem = reshape(collect(1:ref*ref),(ref,ref))
+  slice = transpose(repeat(elem,nC))
+  mat = repeat(slice,nC)
+  return Array(reshape(transpose(mat),nF*nF))
+end
+

src/Refinement/RefinedTriangulations.jl

@@ -0,0 +1,132 @@
+
+
+struct RefinedTriangulation{Dc,Dp,A<:Triangulation,B<:RefinedDiscreteModel} <: Triangulation{Dc,Dp}
+  trian::A
+  model::B
+
+  function RefinedTriangulation(trian::Triangulation{Dc,Dp},model) where {Dc,Dp}
+    A = typeof(trian)
+    B = typeof(model)
+    return new{Dc,Dp,A,B}(trian,model)
+  end
+end
+
+function get_refined_model(t::RefinedTriangulation)
+  return t.model
+end
+
+# Wrap Triangulation API
+function Geometry.get_background_model(t::RefinedTriangulation)
+  get_background_model(t.trian)
+end
+
+function Geometry.get_grid(t::RefinedTriangulation)
+  get_grid(t.trian)
+end
+
+function Geometry.get_glue(t::RefinedTriangulation,::Val{d}) where d
+  get_glue(t.trian,Val(d))
+end
+
+function Base.view(t::RefinedTriangulation,ids::AbstractArray)
+  v = view(t.trian,ids)
+  return RefinedTriangulation(v,t.model)
+end
+
+# Wrap constructors for RefinedDiscreteModel
+function Geometry.Triangulation(
+  ::Type{ReferenceFE{d}},model::RefinedDiscreteModel,filter::AbstractArray) where d
+
+  trian = Triangulation(ReferenceFE{d},get_model(model),filter)
+  return RefinedTriangulation(trian,model)
+end
+
+function Geometry.Triangulation(
+  ::Type{ReferenceFE{d}},model::RefinedDiscreteModel,labels::FaceLabeling;kwargs...) where d
+  trian = Triangulation(ReferenceFE{d},get_model(model),labels;kwargs...)
+  return RefinedTriangulation(trian,model)
+end
+
+function Geometry.Triangulation(trian::RefinedTriangulation,args...;kwargs...)
+  return RefinedTriangulation(Triangulation(trian.trian),trian.model)
+end
+
+# Domain changes
+# TODO: This assumes we have the same type of triangulation on both refinement levels!
+#       we might want to change this in the future when doing hybrid methods etc...
+
+function Geometry.is_change_possible(strian::RefinedTriangulation,ttrian::RefinedTriangulation)
+  (strian === ttrian) && (return true)
+  smodel = get_refined_model(strian)
+  tmodel = get_refined_model(ttrian)
+  a = get_parent(tmodel) === smodel # tmodel = refine(smodel)
+  b = get_parent(smodel) === tmodel # smodel = refine(tmodel)
+  return a || b
+end
+
+function Geometry.is_change_possible(strian::RefinedTriangulation,ttrian::T) where {T <: Triangulation}
+  smodel = get_refined_model(strian)
+  tmodel = get_background_model(ttrian)
+  a = get_model(smodel) === tmodel  # It is fundamentally the same model
+  b = get_parent(smodel) === tmodel # smodel = refine(tmodel)
+  return a || b
+end
+
+function Geometry.is_change_possible(strian::T,ttrian::RefinedTriangulation) where {T <: Triangulation}
+  return is_change_possible(ttrian,strian)
+end
+
+function Geometry.best_target(strian::RefinedTriangulation,ttrian::RefinedTriangulation)
+  @check is_change_possible(strian,ttrian)
+  smodel = get_refined_model(strian)
+  tmodel = get_refined_model(ttrian)
+  a = get_model(smodel) === get_parent(tmodel)
+  b = get_parent(tmodel) === smodel # tmodel = refine(smodel)
+  a || b ? ttrian : strian
+end
+
+function Geometry.best_target(strian::RefinedTriangulation,ttrian::T) where {T <: Triangulation}
+  @check is_change_possible(strian,ttrian)
+  return strian
+end
+
+function Geometry.best_target(strian::T,ttrian::RefinedTriangulation) where {T <: Triangulation}
+  @check is_change_possible(strian,ttrian)
+  return ttrian
+end
+
+
+"""
+  Given a RefinedTriangualtion and a CellField defined on the parent(coarse) mesh, 
+  returns an equivalent CellField on the fine mesh.
+"""
+function change_domain_c2f(f_coarse, ftrian::RefinedTriangulation{Dc,Dp}) where {Dc,Dp}
+  model  = get_refined_model(ftrian)
+  glue   = get_glue(model)
+  if (num_cells(ftrian) != 0)
+    # Coarse field but with fine indexing, i.e 
+    #   f_f2c[i_fine] = f_coarse[coarse_parent(i_fine)]
+    fcell_to_ccell = glue.f2c_faces_map[Dc+1]
+    m = Reindex(get_data(f_coarse))
+    f_f2c = lazy_map(m,fcell_to_ccell)
+
+    # Fine to coarse coordinate map: x_coarse = Φ^(-1)(x_fine)
+    ref_coord_map = get_f2c_ref_coordinate_map(glue)
+
+    # Final map: f_fine(x_fine) = f_f2c ∘ Φ^(-1)(x_fine) = f_coarse(x_coarse)
+    f_fine = lazy_map(∘,f_f2c,ref_coord_map)
+    return GenericCellField(f_fine,ftrian,ReferenceDomain())
+  else
+    f_fine = Fill(Gridap.Fields.ConstantField(0.0),num_cells(ftrian))
+    return GenericCellField(f_fine,ftrian,ReferenceDomain())
+  end
+end
+
+function CellData.change_domain(a::CellField,ttrian::RefinedTriangulation)
+  strian = get_triangulation(a)
+  if strian === ttrian
+    return a
+  end
+  @assert is_change_possible(strian,ttrian)
+  change_domain_c2f(a,ttrian)
+end

src/Refinement/Refinement.jl

@@ -0,0 +1,42 @@
+"""
+
+"""
+module Refinement
+
+using FillArrays
+using LinearAlgebra
+using IterativeSolvers
+using Gridap.Helpers
+using Gridap.Arrays
+using Gridap.Geometry
+using Gridap.FESpaces
+using Gridap.ReferenceFEs
+using Gridap.CellData
+using Gridap.Visualization
+
+import Base: view
+import LinearAlgebra: mul!
+import Gridap.Geometry: get_grid, get_grid_topology, get_face_labeling
+import Gridap.Geometry: Triangulation, is_change_possible, best_target, get_glue, get_background_model
+import Gridap.CellData: change_domain
+
+include("RefinedDiscreteModels.jl")
+include("RefinedTriangulations.jl")
+include("RefinementTransferOperators.jl")
+
+export RefinementGlue
+export get_f2c_ref_cell_map, get_f2c_ref_coordinate_map
+
+export RefinedDiscreteModel
+export get_model, get_parent, get_glue
+export RefinedCartesianDiscreteModel
+
+export RefinedTriangulation
+export Triangulation, is_change_possible, best_target, get_refined_model
+export change_domain_c2f, change_domain
+
+export RefinementTransferOperator, RefinementTransferMap
+export mul!
+
+
+end # module