Quadrature order and faces shapes

[mancellin]

When meshing the boundary with triangles, setting qorder=1 gives me a quadrature with one quadrature point per face.
When meshing the boundary with quadrilaterals, setting qorder=1 gives me 4 points per face.
That is a bit counterintuitive to me. Is it on purpose?

using Inti, Gmsh

const gmsh = Gmsh.gmsh
const mesh_file = "/tmp/generated_mesh.msh"

function generate_cylinder_mesh(radius, h, resolution; triangles=true)
    gmsh.initialize()
    gmsh.model.add("disk_with_thickness")
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", resolution)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", resolution)
    cyl = gmsh.model.occ.addCylinder(0.0, 0.0, 0.0, 0.0, 0.0, h, radius)
    gmsh.model.occ.synchronize()
    if !triangles
        gmsh.model.mesh.setRecombine(2, cyl)
        gmsh.option.setNumber("Mesh.RecombineAll", 1)
    end
    gmsh.model.mesh.generate(2)
    gmsh.write(mesh_file)
    gmsh.finalize()

    msh = Inti.import_mesh(mesh_file)
    Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 3, msh)
    Γ = Inti.boundary(Ω)
    return msh[Γ]
end

triangle_mesh = generate_cylinder_mesh(1.0, 1.0, 0.1, triangles=true)
triangle_quadra = Inti.Quadrature(triangle_mesh; qorder=1)
println(length(triangle_quadra))  # = number of elements in triangle_mesh

quadri_mesh = generate_cylinder_mesh(1.0, 1.0, 0.1, triangles=false)
quadri_quadra = Inti.Quadrature(quadri_mesh; qorder=1)
println(length(quadri_quadra))  # = 4 times the number of elements in quadri_mesh

[tanderson92]

This behavior is explained by _qrule_for_reference_shape() which selects a default quadrature type when the keyword argument qorder is passed absent a specific direction for the choice of quadrature rule to the Quadrature constructor. In the example you provided, the aforementioned routine chooses tensor-product quadrature rules which you may expect to be somewhat less efficient. I do not see any other quadrature rules implemented for a quadrilateral at the moment.

[maltezfaria]

Yes, as @tanderson92 mentioned we use tensor product quadrature rules for quadrilateral domains. But I don't think there is a rule for integrating $\mathbb{P}_{1,\square} = \{1, x, y, xy\}$ which contains only one point (the space we wish to integrate has dimension 4, and one quadrature point only gives you three degrees of freedom). But in general, yes, more efficient (tabulated?) quadrature rules for the square would be welcome 😄

[maltezfaria]

At some point we switched from Fejer quadrature in 1d to the more efficient Gauss-Legendre quadrature in 1d, which fixes this issue in main since now the tensor product quadrature is done using Gauss-Legendre. Thanks for the MWE @mancellin!