Fixed quadrature orders

Project.toml

@@ -32,6 +32,7 @@ BSON = "0.3"
 FileIO = "1"
 FillArrays = "1.13.0"
 Gridap = "0.18"
+GridapDistributed = "0.4.5"
 GridapGmsh = "0.7.2"
 GridapP4est = "0.3.7"
 GridapPETSc = "0.5.0"

src/main.jl

@@ -255,15 +255,13 @@ end
 function _fe_space(::Val{:p},params)
   @notimplementedif space_uses_multigrid(params[:solver])[2]
 
-  k = params[:fespaces][:order_u]
   Ωf = params[:Ωf]
-
-  order = params[:fespaces][:p_order](k)
+  k = params[:fespaces][:p_order]
   space = params[:fespaces][:p_space]
   conformity = params[:fespaces][:p_conformity]
   constraint = params[:fespaces][:p_constraint]
 
-  reffe_p = ReferenceFE(lagrangian,Float64,order;space)
+  reffe_p = ReferenceFE(lagrangian,Float64,k;space)
   if uses_macro_elements(params)
     rrule = params[:fespaces][:rrule]
     reffe_p = Gridap.Adaptivity.MacroReferenceFE(

src/parameters.jl

@@ -308,10 +308,21 @@ 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
+  # Add discretization parameters
   @assert haskey(FLUID_DISCRETIZATIONS[Symbol(poly)],fespaces[:fluid_disc])
   merge!(fespaces,pairs(FLUID_DISCRETIZATIONS[Symbol(poly)][fespaces[:fluid_disc]]))
-
+  fespaces[:p_order] = fespaces[:p_order](fespaces[:order_u])
+
+  # Integration
+  fespaces[:k] = max(fespaces[:order_u],fespaces[:order_j]) # Maximal polynomial degree
+  fespaces[:q] = max( # Quadrature order
+    2*(fespaces[:order_u] - 1),              # UU-Laplacian
+    3*fespaces[:order_u] - 1,                # UU-Convection
+    2*fespaces[:order_j],                    # JJ-Mass
+    fespaces[:order_u] + fespaces[:order_j], # UJ-Coupling
+  )
+  fespaces[:quadratures] = generate_quadratures(poly,fespaces)
+
   return fespaces
 end
 
@@ -364,9 +375,33 @@ function rt_scaling(model,feparams)
   return phi
 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
+    # 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.
+    quads = map(enumerate(fpolys)) do (d,p)
+      if d == D
+        Quadrature(
+          poly,Gridap.Adaptivity.CompositeQuadrature(),feparams[:rrule],qdegree
+        )
+      else
+        Quadrature(p,qdegree)
+      end
+    end
+  else
+    quads = map(p -> Quadrature(p,qdegree), fpolys)
+  end
+
+  return quads
+end
+
 """
 Valid keys for `params[:multigrid]` are the following:
-
 """
 function params_multigrid(params::Dict{Symbol,Any})
   solver = params[:solver]

src/weakforms.jl

@@ -171,9 +171,10 @@ end
 retrieve_fluid_params(params,k) = retrieve_fluid_params(params[:model],params,k)
 
 function retrieve_fluid_params(model,params,k)
-  fluid = params[:fluid]
-  Ωf    = params[:Ωf]
-  dΩf   = Measure(Ωf,2*k)
+  fluid  = params[:fluid]
+  quad3D = params[:fespaces][:quadratures][3]
+  Ωf  = params[:Ωf]
+  dΩf = Measure(Ωf,quad3D)
 
   α, β, γ, σf = fluid[:α], fluid[:β], fluid[:γ], fluid[:σ]
   f = fluid[:f]
@@ -186,10 +187,11 @@ end
 retrieve_solid_params(params,k) = retrieve_solid_params(params[:model],params,k)
 
 function retrieve_solid_params(model,params,k)
-  solid = params[:solid]
+  quad3D = params[:fespaces][:quadratures][3]
+  solid  = params[:solid]
   if solid !== nothing
     Ωs  = params[:Ωs]
-    dΩs = Measure(Ωs,2*k)
+    dΩs = Measure(Ωs,quad3D)
     σs  = solid[:σ]
     return Ωs, dΩs, σs
   else
@@ -200,13 +202,15 @@ end
 retrieve_bcs_params(params,k) = retrieve_bcs_params(params[:model],params,k)
 
 function retrieve_bcs_params(model,params,k)
+  quad3D = params[:fespaces][:quadratures][3]
+  quad2D = params[:fespaces][:quadratures][2]
   bcs = params[:bcs]
 
   params_φ = []
   for i in 1:length(bcs[:φ])
     φ_i = bcs[:φ][i][:value]
     Γ   = _boundary(model,bcs[:φ][i][:domain])
-    dΓ  = Measure(Γ,2*k)
+    dΓ  = Measure(Γ,quad2D)
     n_Γ = get_normal_vector(Γ)
     push!(params_φ,(φ_i,n_Γ,dΓ))
   end
@@ -217,7 +221,7 @@ function retrieve_bcs_params(model,params,k)
     cw_i = bcs[:thin_wall][i][:cw]
     jw_i = bcs[:thin_wall][i][:jw]
     Γ    = _boundary(model,bcs[:thin_wall][i][:domain])
-    dΓ   = Measure(Γ,2*k)
+    dΓ   = Measure(Γ,quad2D)
     n_Γ  = get_normal_vector(Γ)
     push!(params_thin_wall,(τ_i,cw_i,jw_i,n_Γ,dΓ))
   end
@@ -230,22 +234,22 @@ function retrieve_bcs_params(model,params,k)
   for i in 1:length(bcs[:f])
     f_i  = bcs[:f][i][:value]
     Ω_i  = _interior(model,bcs[:f][i][:domain])
-    dΩ_i = Measure(Ω_i,2*k)
+    dΩ_i = Measure(Ω_i,quad3D)
     push!(params_f,(f_i,dΩ_i))
   end
 
   params_B = []
   for i in 1:length(bcs[:B])
     B_i  = bcs[:B][i][:value]
     Ω_i  = _interior(model,bcs[:B][i][:domain])
-    dΩ_i = Measure(Ω_i,2*k)
+    dΩ_i = Measure(Ω_i,quad3D)
     push!(params_f,(γ,B_i,dΩ_i))
   end
 
   params_Λ = []
   for i in 1:length(params[:bcs][:stabilization])
     Λ = _skeleton(model,params[:bcs][:stabilization][i][:domain])
-    dΛ = Measure(Λ,2*k)
+    dΛ = Measure(Λ,quad2D)
     h = get_cell_size(Λ)
     μ = params[:bcs][:stabilization][i][:μ]
     push!(params_Λ,(μ,h,dΛ))

test/seq/hunt_tests.jl

@@ -23,6 +23,16 @@ hunt(
   solver=:julia,
 )
 
+hunt(
+  nc=(10,10),
+  L=1.0,
+  B=(0.,50.,0.),
+  order=3,
+  title="hunt",
+  solver=:julia,
+  simplexify=true,
+)
+
 hunt(
   nc=(4,4),
   np=(1,1),
@@ -63,14 +73,4 @@ hunt(
   kmap_y = 3
 )
 
-hunt(
-  nc=(4,4),
-  L=1.0,
-  B=(0.,50.,0.),
-  order=3,
-  title="hunt",
-  solver=:julia,
-  simplexify=true,
-)
-
 end # module