Many fixes

Project.toml

@@ -8,6 +8,7 @@ BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
@@ -29,6 +30,7 @@ gmsh_jll = "630162c2-fc9b-58b3-9910-8442a8a132e6"
 [compat]
 BSON = "0.3"
 FileIO = "1"
+FillArrays = "1.13.0"
 Gridap = "0.18"
 GridapGmsh = "0.7.2"
 GridapP4est = "0.3.7"

src/Applications/channel.jl

@@ -300,7 +300,7 @@ function channel_mesh(
     model = get_model(mh,1)
   end
   if simplexify
-    model = simplexify(model)
+    model = Gridap.simplexify(model)
   end
   params[:model] = model
   return model

src/Applications/hunt.jl

@@ -296,7 +296,7 @@ function hunt_mesh(
     model = get_model(mh,1)
   end
   if simplexify
-    model = simplexify(model)
+    model = Gridap.simplexify(model)
   end
   params[:model] = model
   return model

src/GridapMHD.jl

@@ -2,9 +2,7 @@ module GridapMHD
 
 using Random
 using LinearAlgebra
-using SparseArrays
-using SparseMatricesCSR
-using BlockArrays
+using SparseArrays, SparseMatricesCSR, BlockArrays, FillArrays
 using ForwardDiff
 
 using FileIO

src/gridap_extras.jl

@@ -1,4 +1,9 @@
 
+# 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)
@@ -43,3 +48,37 @@ function get_mesh_size(m::GridapDistributed.DistributedDiscreteModel)
   return reduce(min,h;init=one(eltype(h)))
 end
 
+# MacroReferenceFEs
+
+function conformity_from_symbol(sym::Symbol)
+  if sym == :H1
+    return H1Conformity()
+  elseif sym == :L2
+    return L2Conformity()
+  else
+    @assert false
+  end
+end
+
+function Gridap.Adaptivity.MacroReferenceFE(
+  rrule::Gridap.Adaptivity.RefinementRule,
+  reffes::AbstractVector{<:Tuple};
+  conformity = Conformity(first(reffes))
+)
+  polys = Gridap.Adaptivity.get_cell_polytopes(rrule)
+  _reffes = map(polys,reffes) do p,r
+    basis, args, kwargs = r
+    ReferenceFE(p,basis,args...;kwargs...)
+  end
+  return Gridap.Adaptivity.MacroReferenceFE(rrule,_reffes;conformity)
+end
+
+function Gridap.Adaptivity.MacroReferenceFE(
+  rrule::Gridap.Adaptivity.RefinementRule,
+  reffe::Tuple;
+  conformity = Conformity(first(reffes))
+)
+  basis, args, kwargs = reffe
+  reffes = ReferenceFE(rrule.ref_grid,basis,args...;kwargs...)
+  return Gridap.Adaptivity.MacroReferenceFE(rrule,reffes;conformity)
+end

src/main.jl

@@ -260,8 +260,9 @@ function _fe_space(::Val{:p},params)
   reffe_p = ReferenceFE(lagrangian,Float64,order;space)
   if uses_macro_elements(params)
     rrule = params[:fespaces][:rrule]
-    subreffes = Fill(reffe_p,Adaptivity.num_subcells(rrule))
-    reffe_p = Gridap.Adaptivity.MacroReferenceFE(rrule,subreffes;conformity)
+    reffe_p = Gridap.Adaptivity.MacroReferenceFE(
+      rrule,reffe_p;conformity = conformity_from_symbol(conformity)
+    )
   end
   params[:fespaces][:reffe_p] = reffe_p
 
@@ -379,9 +380,6 @@ end
 
 # Sub-triangulations
 
-const DiscreteModelTypes = Union{Gridap.DiscreteModel,GridapDistributed.DistributedDiscreteModel}
-const TriangulationTypes = Union{Gridap.Triangulation,GridapDistributed.DistributedTriangulation}
-
 _interior(model,domain::DiscreteModelTypes) = Interior(domain)
 _interior(model,domain::TriangulationTypes) = domain
 _interior(model,domain::Nothing) = Triangulation(model)

src/parameters.jl

@@ -309,45 +309,43 @@ function params_fespaces(params::Dict{Symbol,Any})
   fespaces = _add_optional(params[:fespaces],mandatory,optional,params,"[:fespaces]")
 
   fespaces[:k] = max(fespaces[:order_u],fespaces[:order_j]) # Maximal polynomial degree
-  @assert haskey(FLUID_DISCRETIZATIONS[poly],fespaces[:fluid_disc])
-  merge!(fespaces,FLUID_DISCRETIZATIONS[poly][fespaces[:fluid_disc]])
+  @assert haskey(FLUID_DISCRETIZATIONS[Symbol(poly)],fespaces[:fluid_disc])
+  merge!(fespaces,pairs(FLUID_DISCRETIZATIONS[Symbol(poly)][fespaces[:fluid_disc]]))
 
   return fespaces
 end
 
 """
-  List of supported velocity-pressure pairs for the fluid discretization.
+    const FLUID_DISCRETIZATIONS
+
+List of possible fluid discretizations for the (velocity,pressure) pair.
 """
 const FLUID_DISCRETIZATIONS = (;
-  HEX => FLUID_DISCRETIZATIONS_HEX,
-  TET => FLUID_DISCRETIZATIONS_TET,
-)
-
-const FLUID_DISCRETIZATIONS_HEX = (;
-  :Qk_dPkm1 => (;
-    :p_space => :P,
-    :p_conformity => :L2,
-    :p_order => (k) -> k-1,
-  ),
-  :Qk_Qkm1 => (; # Taylor Hood
-    :p_space => :Q,
-    :p_conformity => :H1,
-    :p_order => (k) -> k-1,
+  :HEX => (; # Hexahedra
+    :Qk_dPkm1 => (;
+      :p_space => :P,
+      :p_conformity => :L2,
+      :p_order => (k) -> k-1,
+    ),
+    :Qk_Qkm1 => (; # Taylor Hood
+      :p_space => :Q,
+      :p_conformity => :H1,
+      :p_order => (k) -> k-1,
+    ),
   ),
-)
-
-const FLUID_DISCRETIZATIONS_TET = (;
-  :Pk_Pkm1 => (; # Taylor Hood
-    :p_space => :P,
-    :p_conformity => :H1,
-    :p_order => (k) -> k-1,
-  ),
-  :SV => (; # Scott-Vogelius (Pk_dPkm1) with macro-elements
-    :p_space => :P,
-    :p_conformity => :L2,
-    :p_order => (k) -> k-1,
-    :rrule => Gridap.Adaptivity.BarycentricRefinementRule(TET)
-  ),        
+  :TET => (; # Tetrahedra
+    :Pk_Pkm1 => (; # Taylor Hood
+      :p_space => :P,
+      :p_conformity => :H1,
+      :p_order => (k) -> k-1,
+    ),
+    :SV => (; # Scott-Vogelius (Pk_dPkm1) with macro-elements
+      :p_space => :P,
+      :p_conformity => :L2,
+      :p_order => (k) -> k-1,
+      :rrule => Gridap.Adaptivity.BarycentricRefinementRule(TET)
+    ),        
+  )
 )
 
 function uses_macro_elements(params::Dict)
@@ -535,7 +533,7 @@ function params_bcs_u(params::Dict{Symbol,Any})
   )
   _check_mandatory(params[:bcs][:u],mandatory,"[:bcs][:u]")
   optional = Dict(
-    :values=>zero_values(params[:bcs][:u][:tags]),
+    :values => zero_values(params[:bcs][:u][:tags]),
   )
   u = _add_optional(params[:bcs][:u],mandatory,optional,params,"[:bcs][:u]")
   _check_unused(u,mandatory,params,"[:bcs][:u]")

test/seq/hunt_tests.jl

@@ -63,4 +63,14 @@ hunt(
   kmap_y = 3
 )
 
+hunt(
+  nc=(2,2),
+  L=1.0,
+  B=(0.,50.,0.),
+  order=3,
+  title="hunt",
+  solver=:julia,
+  simplexify=true,
+)
+
 end # module