Minor refactor to accomodate HDiv-conforming velocity

src/Applications/channel.jl

@@ -160,7 +160,7 @@ function _channel(;
 
   # Fluid parameters
   params[:fluid] = Dict(
-    :domain=>nothing,
+    :domain => nothing,
     :α=>α,
     :β=>β,
     :γ=>γ,

src/Applications/hunt.jl

@@ -287,9 +287,9 @@ function hunt_mesh(
   ranks_per_level = nothing, adaptivity_method = 0
 )
   if isnothing(ranks_per_level) # Single grid
-    model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted,mesh_postpro) 
+    model = Meshers.hunt_generate_base_mesh(parts,np,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted) 
   else # Multigrid
-    mh = Meshers.hunt_generate_mesh_hierarchy(parts,ranks_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted,mesh_postpro)
+    mh = Meshers.hunt_generate_mesh_hierarchy(parts,ranks_per_level,nc,L,tw,Ha,kmap_x,kmap_y,BL_adapted)
     params[:multigrid] = Dict{Symbol,Any}(
       :mh => mh,
       :num_refs_coarse => 0,

src/GridapMHD.jl

@@ -40,6 +40,7 @@ include("Solvers/badia2024.jl")
 # Main driver
 include("gridap_extras.jl")
 include("utils.jl")
+include("geometry.jl")
 include("parameters.jl")
 include("weakforms.jl")
 include("main.jl")

src/Solvers/li2019.jl

@@ -16,7 +16,7 @@ function Li2019Solver(op::FEOperator,params)
   # Preconditioner
   model = params[:model]
   k  = max(params[:fespaces][:order_u],params[:fespaces][:order_j])
-  Ωf = _interior(model,params[:fluid][:domain])
+  Ωf = params[:Ωf]
   dΩ = Measure(Ωf,2*k)
   a_Dj = li2019_Dj(dΩ,params)
   a_Fk = li2019_Fk(dΩ,params)

src/geometry.jl

@@ -0,0 +1,148 @@
+
+const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel}
+const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation}
+
+function setup_geometry!(params)
+  params[:geometry] = Dict{Symbol,Any}(
+    :interiors  => Dict{UInt64,Any}(),
+    :boundaries => Dict{UInt64,Any}(),
+    :skeletons  => Dict{UInt64,Any}(),
+    :normals    => Dict{UInt64,Any}(),
+    :measures   => Dict{UInt64,Any}(),
+    :other      => Dict{UInt64,Any}(),
+  )
+
+  Ω = interior(params,nothing)
+  params[:Ω] = Ω
+
+  fluid = params[:fluid]
+  Ωf = interior(params,fluid[:domain])
+  params[:Ωf] = Ωf
+
+  solid = params[:solid]
+  Ωs = !isnothing(solid) ? interior(params,solid[:domain]) : nothing
+  params[:Ωs] = Ωs
+
+  if uses_multigrid(params[:solver])
+    params[:multigrid][:Ω]  = params[:multigrid][:mh]
+    params[:multigrid][:Ωf] = params[:multigrid][:mh]
+    params[:multigrid][:Ωs] = params[:multigrid][:mh]
+  end
+end
+
+domain_hash(domain::Union{Nothing,DiscreteModelTypes,TriangulationTypes}) = objectid(domain)
+domain_hash(domain::Union{String,Vector{String}}) = hash(join(domain))
+
+domain_hash(model,domain::Union{String,Vector{String}}) = hash(join(domain),objectid(model))
+domain_hash(model,domain) = domain_hash(domain)
+
+interior(params::Dict,domain) = interior(params,params[:model],domain)
+
+function interior(params::Dict,model,domain)
+  interiors = params[:geometry][:interiors]
+  key = domain_hash(model,domain)
+  if haskey(interiors,key)
+    trian = interiors[key]
+  else
+    trian = _interior(model,domain)
+    interiors[domain_hash(trian)] = trian
+    interiors[key] = trian
+  end
+  return trian
+end
+
+boundary(params::Dict,domain) = boundary(params,params[:model],domain)
+
+function boundary(params::Dict,model,domain)
+  boundaries = params[:geometry][:boundaries]
+  key = domain_hash(model,domain)
+  if haskey(boundaries,key)
+    trian = boundaries[key]
+  else
+    trian = _boundary(model,domain)
+    boundaries[domain_hash(trian)] = trian
+    boundaries[key] = trian
+  end
+  return trian
+end
+
+skeleton(params::Dict,domain) = skeleton(params,params[:model],domain)
+
+function skeleton(params::Dict,model,domain)
+  skeletons = params[:geometry][:skeletons]
+  key = domain_hash(model,domain)
+  if haskey(skeletons,key)
+    trian = skeletons[key]
+  else
+    trian = _skeleton(model,domain)
+    skeletons[domain_hash(trian)] = trian
+    skeletons[key] = trian
+  end
+  return trian
+end
+
+function normal_vector(params::Dict,trian::TriangulationTypes)
+  normals = params[:geometry][:normals]
+  key = domain_hash(trian)
+  if haskey(normals,key)
+    n = normals[key]
+  else
+    n = get_normal_vector(trian)
+    normals[key] = n
+  end
+  return n
+end
+
+function measure(params::Dict,trian::TriangulationTypes)
+  measures = params[:geometry][:measures]
+  key = domain_hash(trian)
+  if haskey(measures,key)
+    m = measures[key]
+  else
+    d = num_cell_dims(trian)
+    q = params[:fespaces][:quadratures][d]
+    m = Measure(trian,q)
+    measures[key] = m
+  end
+  return m
+end
+
+_interior(model,domain::DiscreteModelTypes) = Triangulation(domain)
+_interior(model,domain::TriangulationTypes) = domain
+_interior(model,domain::Nothing) = Triangulation(model)
+_interior(model,domain) = Triangulation(model,tags=domain)
+
+_boundary(model,domain::TriangulationTypes) = domain
+_boundary(model,domain::Nothing) = Boundary(model)
+_boundary(model,domain) = Boundary(model,tags=domain)
+
+_skeleton(model,domain::TriangulationTypes) = Skeleton(domain)
+_skeleton(model,domain::Nothing) = Skeleton(model)
+_skeleton(model,domain) = Skeleton(_interior(model,domain))
+
+# Face/cell sizes
+
+get_cell_size(t::TriangulationTypes) = CellField(_get_cell_size(t),t)
+
+get_cell_size(m::DiscreteModelTypes) = get_cell_size(Triangulation(m))
+_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m))
+
+function _get_cell_size(t::Triangulation) :: Vector{Float64}
+  if iszero(num_cells(t))
+    return Float64[]
+  end
+  meas = get_cell_measure(t)
+  d = num_dims(t)
+  return collect(Float64, meas .^ (1/d))
+end
+
+function _get_cell_size(t::GridapDistributed.DistributedTriangulation)
+  map(_get_cell_size,local_views(t))
+end
+
+get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m))
+
+function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel)
+  h = map(get_mesh_size,local_views(m))
+  return reduce(min,h;init=one(eltype(h)))
+end

src/gridap_extras.jl

@@ -1,9 +1,4 @@
 
-# Quality of life definitions
-
-const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel}
-const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation}
-
 # Get polytopes
 
 function Geometry.get_polytopes(model::GridapDistributed.DistributedDiscreteModel)
@@ -21,33 +16,6 @@ function Geometry.get_polytopes(trian::GridapDistributed.DistributedTriangulatio
   return getany(polys)
 end
 
-# Face/cell sizes
-
-get_cell_size(t::TriangulationTypes) = CellField(_get_cell_size(t),t)
-
-get_cell_size(m::DiscreteModelTypes) = get_cell_size(Triangulation(m))
-_get_cell_size(m::DiscreteModelTypes) = _get_cell_size(Triangulation(m))
-
-function _get_cell_size(t::Triangulation) :: Vector{Float64}
-  if iszero(num_cells(t))
-    return Float64[]
-  end
-  meas = get_cell_measure(t)
-  d = num_dims(t)
-  return collect(Float64, meas .^ (1/d))
-end
-
-function _get_cell_size(t::GridapDistributed.DistributedTriangulation)
-  map(_get_cell_size,local_views(t))
-end
-
-get_mesh_size(m::DiscreteModel) = minimum(_get_cell_size(m))
-
-function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel)
-  h = map(get_mesh_size,local_views(m))
-  return reduce(min,h;init=one(eltype(h)))
-end
-
 # MacroReferenceFEs
 
 function conformity_from_symbol(sym::Symbol)
@@ -91,6 +59,3 @@ function Gridap.Adaptivity.MacroReferenceFE(
   reffes = Fill(reffe,Gridap.Adaptivity.num_subcells(rrule))
   return Gridap.Adaptivity.MacroReferenceFE(rrule,reffes;macro_kwargs...)
 end
-
-# Local projection operators
-

src/main.jl

@@ -124,8 +124,7 @@ function main(_params::Dict;output::Dict=Dict{Symbol,Any}())
   params = add_default_params(_params)
   t = params[:ptimer]
 
-  # Compute triangulations only once for performance
-  _setup_trians!(params)
+  setup_geometry!(params)
 
   # FESpaces
   tic!(t;barrier=true)
@@ -256,7 +255,6 @@ function _fe_space(::Val{:p},params)
 end
 
 function _fe_space(::Val{:j},params)
-  k = params[:fespaces][:order_j]
   uses_mg = space_uses_multigrid(params[:solver])[3]
   model = uses_mg ? params[:multigrid][:mh] : params[:model]
 
@@ -355,51 +353,6 @@ function _continuation_fe_operator(mfs,U,V,params)
   return op
 end
 
-# Sub-triangulations
-
-_interior(model,domain::DiscreteModelTypes) = Interior(domain)
-_interior(model,domain::TriangulationTypes) = domain
-_interior(model,domain::Nothing) = Triangulation(model)
-_interior(model,domain) = Interior(model,tags=domain)
-
-_boundary(model,domain::TriangulationTypes) = domain
-_boundary(model,domain::Nothing) = Boundary(model)
-_boundary(model,domain) = Boundary(model,tags=domain)
-
-_skeleton(model,domain::TriangulationTypes) = SkeletonTriangulation(domain)
-_skeleton(model,domain::Nothing) = SkeletonTriangulation(model)
-_skeleton(model,domain) = _skeleton(model,_interior(model,domain))
-
-function _setup_trians!(params)
-  Ω = _interior(params[:model],nothing)
-
-  fluid = params[:fluid]
-  Ωf = _interior(params[:model],fluid[:domain])
-
-  solid = params[:solid]
-  Ωs = !isnothing(solid) ? _interior(params[:model],solid[:domain]) : nothing
-
-  #if !uses_multigrid(params[:solver])
-    params[:Ω]  = Ω
-    params[:Ωf] = Ωf
-    params[:Ωs] = Ωs
-  #else
-  #  params[:multigrid][:Ω]  = Ω
-  #  params[:multigrid][:Ωf] = Ωf
-  #  params[:multigrid][:Ωs] = Ωs
-  #  params[:Ω]  = Ω[1]
-  #  params[:Ωf] = Ωf[1]
-  #  params[:Ωs] = Ωs[1]
-  #end
-
-  if uses_multigrid(params[:solver])
-    params[:multigrid][:Ω]  = params[:multigrid][:mh]
-    params[:multigrid][:Ωf] = params[:multigrid][:mh]
-    params[:multigrid][:Ωs] = params[:multigrid][:mh]
-  end
-
-end
-
 # Random vector generation
 
 function _rand(vt::Type{<:Vector{T}},r::AbstractUnitRange) where T

src/parameters.jl

@@ -372,6 +372,8 @@ const FLUID_DISCRETIZATIONS = (;
   )
 )
 
+has_hdiv_fluid_disc(params) = params[:fespaces][:fluid_disc] ∈ (:RT, :BDM)
+
 function fluid_discretization(disc::Symbol,poly::Polytope,feparams)
   D = num_dims(poly)
   k = feparams[:order_u]
@@ -509,11 +511,10 @@ function rt_scaling(poly,feparams)
 end
 
 function generate_quadratures(poly::Polytope{D},feparams) where D
-  fluid_disc = feparams[:fluid_disc]
   qdegree = feparams[:q]
 
   fpolys = ReferenceFEs.get_reffaces(poly)[2:end] # Skip 0-dfaces
-  if fluid_disc == :SV
+  if haskey(feparams,:rrule)
     # TODO: This should be made more general, to ensure all d-faces are properly 
     # refined if need be. In the case of SV, there is not subdivision of the 
     # facets and edges, so its fine.